Nonergodicity parameters for a molecular liquid: 
A comparison between mode coupling theory and simulation 
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We apply the mode coupling theory (MCT) which was recently worked out for molecular liquids 
to a liquid of diatomic, rigid molecules. Making use of the static correlators S™,(q) from a molec- 
ular dynamical simulation, we have solved the MCT-equations for the nonergodicity parameters 
fi n (q) = linit^oo S™(q,t)/S™(q), assuming that all correlators are nonzero for I = I' = 0, 1, ...,/ co , 
only. Depending on l co = 0, 1,2 we obtain different types of ideal glass transitions with transition 
temperatures T c which are below T^ 10 as deduced from the MD-simulation. The q-dependence of 
the critical nonergodicity parameter {f™)^? from the simulation is reasonably well reproduced by 
the corresponding MCT-result. The influence of a strong and weak diagonalization approximation 
on the results is discussed. 
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I. INTRODUCTION 

The use of mode coupling theory (MCT) is one of the 
most important attempts to describe the structural glass 
transition. Restricted to simple liquids, it has provided 
an equation of motion for the normalized density corre- 
lator 



*(<?,*) 



S(q,t) 
S(q) 



(1) 



where S(q) is the corresponding static correlation func- 
tion. MCT, in its idealized version, interprets the glass 
transition as a dynamical transition at a critical temper- 
ature T c (or critical density n c ) from a ergodic (super- 
cooled liquid) to a nonergodic phase (glass). In the lat- 
ter phase the density p(q, t) does not approach its equi- 
librium value po = (p(q,t)) m the long time limit, i.e. 
the correlation function <3>(g, t) of the density fluctuations 
Sp(q, t) — p(q, t) — po does not decay to zero for t — > oo. 
Therefor one can use the nonergodicity parameter (NEP) 



f(q) = lim <f>(q,t) 



(2) 



as an order parameter for the glass transition. f{q) van- 
ishes above and becomes nonzero at T c . In the long time 
limit t ^> to the MCT-equations involve only S(q) as an 
input quantity. Here, to is a microscopic time of order 
10~ 13 sec. This fact demonstrates the microscopic na- 
ture of MCT, because S(q) can be calculated from the 
potential energy V(xi, x^v) which only depends on the 
microscopic configuration { x„ } of the particles. The 
discussion of the MCT-equations for T in the vicinity 
of T c leads to several interesting predictions for the time 
and temperature dependence of <&(<?,£). However, due 
to the existence of crgodicity restoring processes in real 
glass formers below T c , which are taken into account by 
the extended MCT, the transition is smeared out. Never- 
theless, a couple of experimental and numerical investi- 



gations have clearly shown the validity of several predic- 
tions of idealized MCT. For more details the reader may 
consult the reviews [0-1). 

These tests of idealized MCT (simple liquids) were 
mainly qualitative, i.e. the validity of the asymptotic 
power laws for the a- and /3-relaxation regimes was stud- 
ied. Only very few quantitative test exist, which is, be- 
sides the calculation of the critical temperature (or criti- 
cal density), the q-dependence of the critical NEP f c (q) 
at T c (or n c ) and the critical amplitude h(q). This was 
done for a monoatomic (t]] (only T c was determined) and 
binary Lennard- Jones liquid for a hard sphere |J 
and a soft sphere system |^],[l0| and for the molecular liq- 
uids orthoterphenyl jll| and water |jl2|| . Although the 
two latter liquids contain orientational degrees of free- 
dom (ODOF), the analysis of the experimental jll| and 
numerical results [Q was done in the framework of MCT 
for simple liquids, which does not account for the ODOF. 
A first attempt to include ODOF and their coupling to 
the translational degrees of freedom (TDOF) has recently 
been made for water |0| . These authors suggest a semi- 
schematic model where the full q-dependence is consid- 
ered, but the influence of the ODOF is condensed in a 
single parameter \R- The phcnomcnological introduc- 
tion of xr can be viewed as a "renormalisation" of the 
mode coupling strength for the TDOF (center of mass). 

Because most of the glass formers are of molecu- 
lar origin, it is important to extend MCT to molecu- 
lar liquids. This has recently been done for a single 
dumb-bell molecule in an isotropic liqui d |l4|], for linear 
molecules (jl5| and arbitrary molecules |16|-|18[. In the 
present paper we will perform a quantitative test of the 
molecular MCT for a liquid of rigid diatomic molecules 
by comparing the corresponding nonergodicity parame- 
ters as obtained from a recent MD-simulation with those 
from molecular MCT. 

The outline of our paper is as follows. The next section 
contains a descripition of the molecular MCT-equations 



1 



for linear molecules. The diatomic liquid and the com- 
putational details are described in section 3. Section 4 
contains the comparison between MCT and numerical re- 
sults and the final section presents a summary and some 
conclusions. 



with 



S(q,z) = - zl + k(q,z)S~\q) S(q) (7) 



II. MOLECULAR MCT 

In this section we give a short summary of the ap- 
plication of MCT to molecular liquids. We will restrict 
ourselves to linear molecules. For details the reader is 
referred to Ref. The starting point is the tensorial 

one-particle density 

N 

Pim(q, t) = V4tt i l ^2 exp(iq x n (t)) Yi m (tt n (t)) , 
I = 0,1,2,... ,-l<m<l (3) 



where x n (t) and Q n (t) = (9 n (t),4> n (t)) are, respectively, 
the center of mass position and the Euler angles of the 
n-th molecule at time t. Due to the axial symmetry of 
the molecule the third Euler angle Xn{t) does not play 
a role. N is the number of molecules and Yi m are the 
spherical harmonics. The factor in front of the sum is 
for technical convenience. The generalization of the den- 
sity correlator S(q 7 t) (cf. eq.(|l])) to molecular liquids of 
linear molecules is straightforward: 



Si m ,i'm'(q,t) = (p* m (q,t)pi> m >(q,0)) 



(4) 



where * denotes complex conjugation and the angular 
brackets refer to the canonical average over the initial 
conditions in phase space. Generally, the correlation ma- 
trix S(q, t) = (Si mt i' m i(q, t)) is nondiagonal in both I and 
m. In the q-frame, i.e. for q = q Q = (0,0,9), Q = |?|) 
S(q, t) becomes diagonal in m: 



II 



(5) 



The matrix elements Sf^q^t) are real, and nonzero for 
|m| < min(7,r), only. In addition they fulfill 15 



< 



S„r(q,t) = S%(q,t). 



(6) 



In the third and fourth section we will restrict ourselves 
to those correlators with I = V = 0,1 and 2. In that case, 
only six independent correlators S" l (q,t) = S™(q,t), 
< m < I exist. 

A closed set of equations can be derived for the Laplace 
transform S(q, z) of S(q,i)by use of the Zwanzig-Mori 
projection formalism in combination with the mode cou- 
pling approximation which replaces the slow part of the 
time dependent memory kernel by a linear combination 
of bilinear terms S iimi ^ m ^ (q 1; t)Si 27n2y i> 27n > 2 (q 2 ,t). The 
final result for the MCT-equations is as follows |D|: 



and 



q?(q) 



k(q,z) = 



a =T 
a =R 



zS- 1 



J- 1 M( 9 ,z)J- 



(9) 



(10) 



= ( T aa> 1 



is the correlation matrix of the transla- 



tional (a — T) and rotational (a = R) current density 
jfrn(l) El- The kernel J _1 M(q, t)J _1 cannot be calcu- 
lated exactly. The mode coupling approximation yields 
for its slow part (which determines the long time relax- 
ation): 



(J^M^J- 
1 

= 2N 



* AV J slow 
E ^ a '(gAA'|g 1 A 1 Ai;q 2 A 2 A^)x 

xS XlK ( qi ,t)S X2X2 (q 2 ,t). (11) 



Here we have used the short hand notation A = (l,m). 
The vertices V aa only depend on the static correlators 
S\\i (q) or cquivalently on the direct correlations func- 
tions cw(q). The explicit expression for V aa can be 
found in Ref. Q. Eqs.(|7|)-(|ll]) are the closed set of 
MCT-equations for linear molecules, and their structure 
is quite similar to that of the MCT-equations for simple 
liquids. But due to the presence of ODOF and their cou- 
pling to TDOF, which is in the vertex functions V aa , the 
mode coupling polynomial ([ll]) is much more involved. 
Therefore, in a first step we will simplify these equations 
by assuming diagonality of S(q,t) and m(q, t) with re- 
spect to I. We stress that this diagonalization approxi- 
mation is also assumed for the static correlators, in order 
to keep the coupled set of equations as simple as possible. 
Since J is already diagonal, eqs.(^)-(^T]) allow for a self 
consistent solution for the normalized correlators 



SIT (3) 



(12) 



in the q-frame. Then, the corresponding nonergodicity 
parameters 

fr(q) = lim SJ»(g,i) = - lim z^^z) (13) 
are solutions of the following set of coupled equations: 
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with the 2x2 matrix 



{F{n) TT {Fr) TH 



In the thermodynamic limit it is: 
(FP(q;{fP(q)})) aa ' = 

d qi / dq 2 (H?( q ,q u q 2 ;{fr(q)})) c 

where 



will exhibit a type-B transition and those with odd I 
may show a type-A or type-B transition. In case of a 
type-A transition the corresponding transition tempera- 

(14) ture is below that of the type-B transition for f{ n with I 
even. Second, whereas eq.(|l4|) for /° contains on its r.h.s. 
/o% > A mi /r, f2 ni fP 2 etc., the equation for, e.g. f? 
involves /rVr, /f 1 fofP etc., but not / %°, due 

(15) to the diagonalization approximation of the static corre- 
lators. Therefore the freezing of the center of mass order 
parameter /q does not necessarily imply the freezing of 
/™. /™ may freeze at a lower temperature via a type-B 
transition. Without the diagonalization approximation 
with respect to the Vs, we expect all fp(q) to freeze via 

(16) a type-B transition at a single critical temperature T c . 
For more details on type-A and type-B transitions the 
reader is referred to the schematic models which are dis- 
cussed in ref. ifO] and in ref . |19f| . 



(Hr(<?,gi, 92 ;{/r(<z)})r = 

= £' E y aa '(l lm \qihm 1 ;q 2 l 2 m 2 )f^(q 1 )f^(q 2 ). (17) 



lxl 2 mini] 



l denotes the summation over l\ and l 2 such that 



h + l 2 + I is even 
\h-h\ <l<h + h- 



(18) 



This restriction results from the Clebsch-Gordon coeffi- 
cients C(lil 2 l;000) which enter into v aa . The expres- 
sions for the v aa are rather involved. Therefore we do 
not give them explicitely, here. They can be deduced 
from V aa (see second paper of ref. by use of: 

i) transformation to the q-frame where one has to use the 
transformation properties of the correlators under rota- 
tions (see again second paper of ref. Jl5|) 

ii) the diagonalization approximation with respect to 1 
and 

iii) the thermodynamic limit N — > oo , V — > oo such that 
N/V = const, which replaces the sums over <7i,Q2m 
eq. (|Tl]) by integrals; then taking into account the con- 
servation of momentum q = q± + q-2 and performing the 
angular integration one obtains the integral representa- 
tion of cq.(|16|) 

yaa _ This type of approximation will be called weak 
diagonalization approximation (WDA) . Assuming in ad- 
dition the diagonality of v aa with respect to a, and ac- 
cordingly the diagonality of F™ (which we call the strong 
diagonalization approximation (SDA)), the eqs.(^)-(|l7|) 
simplify even more. It is the SDA which has been used 
for the calculation of //"(?) for a dumb-bell molecule in 
a liquid of hard spheres jl4| and for a system of dipolar 
hard spheres ]T^ ]. The diagonalization approximation 
for the time-dependent and the static correlators with 
respect to the Vs implies (JlS ) . This restriction has two 
consequences. First, the eq.(14|) for the /" l with / even 
may contain on its r.h.s. bilinear terms /; m // m , but for 
I odd this cannot happen. Therefore /™ with I even 



III. MODEL 



In order to test 
MD-simulation for a 
molecules 



20-22 



the MCT-predictions, a 
system of diatomic, rigid 
was performed. Each molecule is com- 
posed of two different Lennard- Jones particles A and 
B, which are separated by a distance d. Both parti- 
cles have the same mass m. The interaction between two 
molecules is given by the sum of the interactions between 
the four particles which are given by the Lennard- Jones 
potentials V a p{r) = 4e Q/3 [{o- a p/r) 12 - {a aj3 /rf] . The 
Lennard-Jones parameters are given by ctaa = &ab = 1 , 
o~bb = 0.95, £aa — £ab — 1 and ebb = 0.8. We use 
lengths in units of a aa and is taken as unit for the 
temperature (ks — 1). In these units we have found that 
d = 0.5 is an appropriate value in order to avoid forma- 
tion of a liquid crystalline phase. The MD-simulation 
was done for N = 500 molecules in an (N,p,T)-ensemble 
with p — 1.0 and T between 0.477 and 5.0 by use of the 
rattle algorithm. The length of the runs was long enough 
in order to prepare a well equilibrated system. In order 
to increase the statistics, average over at least eight inde- 
pendent runs was performed. One of the main results we 
have found is the existence of a type-B transition with a 
single transition temperature T^ ID = 0.475, at which the 
dynamics crosses over from ergodic to quasi-nonergodic 
behavior. This temperature has been deduced from the 
a-relaxation time of the purely orientational (i.e. 5 = 
in eq. fj|) EQ ], purely translational (i.e. / = /' = in 
eq.P)) |2l|] and of the general correlators where g ^ 
and ( =7 ^ p2| . The critical temperature T c can 
be located within about 2% from these data. For more 
details the reader is referred to refs. [Eo| |22l. 
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IV. RESULTS 

As already mentioned in section |l| we need for the so- 
lution of eqs.(|l4|)-(|l7|) the static correlators S%(q) as a 
function of (q, I, V , m) and temperature T. We remind the 
reader that we use the q-frame representation, without 
restricting generality. The MD-rcsult for S™, (q) is shown 
in Figure 1 for T = 0.49 and for 1 = 1' = 0,1, 2. The vari- 



ation of Sj n (q) = S™(q) with temperature is smooth. 
From this Figure we observe that the q-dependence of 
SJ n (q) strongly depends on I and m. Whereas S^q) 
and Si(q) resembles the behavior of S(q) for a simple 
liquid, with e.g. a well pronounced main maximum at 
qmax ~ 27T, the q-dependence of the other correlators is 
qualitatively different. For instance S% has its main peak 
at q = 0. 
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FIG. 1. ST(q) versus q for T = 0.49: (a) I = m = 0, (b) / = 1, m = 0, 1 and (c) I = 2,m = 0, 1,2. 



As will be seen below, we need SJ n (q) for temper- 
atures below T min = 0.477, the lowest value of our 
MD-simulation. This has been achieved by a linear ex- 
trapolation from the higher temperature regime (T > 
0.477). After smoothing the q-dependence of SJ n (q) 
by a spline under tension, equations (O)— (17) were 
solved by iteration. This procedure is standard and 
had been used for simple liquids (see e.g. pj and M). 
The iteration was stopped after the maximum value of 

\(fr(i)y u+1) - (fri^y^l was less than io ~ 6 - We have 

used SDA and WDA (see section 0). The /-dependence 
has been truncated at l co = 0,1 and 2. The integrations 



in ( |l6| ) were performed numerically. For this the interval 
[0, q co ] was divided into 300 equidistant grid points. Since 
the MD-data were limited to q < q co = 25, we were not 
able to check the sensitivity of the results for /™ (q) on 
a variation with q co . But there should not be a consid- 
erable change by extending q co to higher values than 25, 
because S] n (q) already approaches one for q = 25 (for 
S22 this is not quite true) (cf. Fig.l). 

Special attention must be paid to g = 0. Since 
(Fl n (q,{r i n {q)})) aa ' for I ^ is proportional to q~ 2 
and g _1 for a = a' — T and a 7^ a', respectively, 
we have treated eq. dT3) separately for q — and 
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q > Imin — 0.08. The results for / z m (g) in between were 
obtained by a polynomial interpolation of degree four. 

The first observation we make is that = for 

T > 0.477. This result, which holds for SDA, WDA and 
l co = 0, 1,2, is unexpected, because for all known cases 
MCT overestimates the transition temperature T c . For 
instance, it is T C MCT « 2T C MD for the binary Lcnnard- 
Jones liquid §. For T < T C MI? = 0.475 ideal glass tran- 
sitions are found for our molecular liquid. The transition 
temperatures are given in Table 1. Here a comment is 
in order. The transition temperatures are determined 
whitin an accuracy of 10~ 3 . Since the critical NEP vary 
as (T — T c )2 they are correct only within 3-5 %. For an 
accuracy of 1 % and less one must determine T c more pre- 
cisely. This can be achieved by use of the T-dependencc 
of the largest eigenvalue Eq (T) of the stability matrix |Q , 
since (1 - E {T)) 2 ~ T c - T for T < T c . 

From this table we recognize that the center of mass 
modes (I = m = 0) freeze at T l c =0 = 0.383 independent 
on l co and SDA or WDA. This transition is of type B. 
Ideal glass transitions for I — 1 and 2 occur at lower tem- 
peratures. Their transition scenario depends on whether 
lco = 1 or l co — 2. In case that l co = 1, i.e. all / > 2 are 
neglected, a type-A transition of /{"(<?) exists at 0.192 
for WDA and at a slightly lower value of 0.187 for SDA. 
But for l co = 2 both modes with / = 1 and I = 2 simulta- 
neously freeze at T c z=1 = T l = 2 = 0.310(0.256) for WDA 
(SDA). This transition is of type B for / = 2 and 1 = 1. 
Whereas I = 2 must freeze via a discontinuous (type-B) 



transition, I = 1 could undergo a continuous transition 
(type-A). The coupling between the modes with I = 1 
and I = 2 in eq. (|lj]) is obviously large enough to lead to a 
type-B transition for I = 1, too. These possible scenarios 
are in accordance with the discussion at the end of section 
H. Since for WDA the r.h.s. of eq.([l4])contains more con- 
tributions from the coupling between ODOF and TDOF 
than for SDA, it is reasonable that the critical tempera- 
tures for WDA are larger (or equal) than those for SDA. 

The NEP (//")f D obtained from the MD- 
simulation p2[ by a fit with the von Schweidler law 
and (/ i m ) M * 7 from molecular MCT for l co = 2, WDA 
and T = Ti =1 = T l = 2 = 0.310 are shown in Figure 2. 



Whereas (fj 71 )^ for I 



0,1, 2 and (/r)f CT for/ 

.MCT 



1,2 

are the critical NEP, this is not true for (/q J , since 
Tj: =a > 0.310. Despite the approximations we made, the 
q-dependence of (fj 71 )^ 10 is reasonably well reproduced 

by the MCT-result. Even the shoulder in (fo)^ ID at 
q w 3 is reproduced as a small hump at that q-value. 
Apart from this prepeak, the q-variation of f™ is in 
phase with that of the corresponding static correlators 
(cf. Fig.l). Whereas (fo) MCT > (fo)^ D (as expected 
again, due to the mean field character of MCT) the oppo- 
site is true for / ; m with I = 1 and 2. For SDA the result 

\MCT 



at T = Ti =1 = T, 



for {ff 
in Figure 2) is about 5 
result for WDA. 



= 0.256 (not shown 
10 % below the corresponding 





rpl=0 

■*■ c 


rpl=\ 

c 




lco = 


0.383 
(0.383) 






lco — 1 


0.383 
(0.383) 


0.192 
(0.187) 




/ — 9 

<<co — & 


0.383 
(0.383) 


0.310 
(0.256) 


0.310 
(0.256) 



TABLE I. Critical temperatures T l c for l co = 0, 1, 2 and WDA. The corresponding values for SDA are given in parantheses. 
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FIG. 3. Comparison of (f§) 



(/o U 



0\MCT rrpl=0 



{Tt° = 0.383) and (fo) (T 0.310) 



Finally we compare (/^)f CT at T = T l =° = 0.383, 
i.e. the critical NEP, with that at T w 0.310. The re- 
sult, including (fg) is given in Figure 3. Now, we 
also find that (/°)f ° T < (f§)? D . That (/°)f OT < 

(/o) ( T = 0.310) is trivial, because /, m increases 
monotonously with decreasing temperature. 



V. SUMMARY AND CONCLUSIONS 

For the first time, we have compared numerical results 
for a molecular liquid of diatomic, rigid molecules 20 2^] 
with the corresponding results from molecular MCT for 
linear molecules, which was worked out recently |]l5| , |l8| . 
This comparison was done for the nonergodicity param- 
eters /f(g). First of all, we find that MCT does not 
yield an ideal glass transition in the temperature range 
of the simulation. This is quite similar to what has been 
found for the MD-simulation of water JlJ,[l3| , but is in 
variance with the behavior for simple Lennard-Jones liq- 
uids |7],|8j and hard spheres [QJ|] . For the hard sphere sys- 
tem this means that the critical density n c from MCT is 
less that the experimental value. This difference between 
simple and molecular liquids probably reveals the role of 
ODOF and their coupling to TDOF. This coupling is un- 
derestimated by both diagonalization approximations we 
used. With these approximations the ideal glass transi- 
tion occurs below T^ 1 , only. The existence of two dif- 
ferent transition temperatures for l co — 1 and l co = 2 we 
believe is an artifact of the diagonalization approxima- 
tion (with respect to I). Nevertheless, the q-dependence 
of (/; m )* /£> is reasonably well reproduced by the corre- 
sponding MCT-result. Whether the hump in (/o)* /CT 

at q w 3, which relates to the shoulder of (fo)^, ID at 
the same q-value is genuine can not presently be de- 
cided. But it is interesting that a similar hump in the 
critical NEP has recently been found for the molecular 



glass former orthoterphenyl p3| . Whether its existence 
in the nonergodicity parameter (f§(q)) for the center of 
mass motion is an indirect influence of the ODOF or not 
remains an open question. For OTP this prepeak may 
be related to an intramolecular property, not related to 
ODOF j2§. 

To obtain the MCT-results presented in this paper al- 
ready a considerable numerical effort was necessary. Nev- 
ertheless it will be important to solve the MCT-equations 
without any diagonalization approximation, which we be- 
lieve guarantees (for our molecules) the existence of a sin- 
gle transition temperature T^ CT . In this respect we also 
would like to mention that, e.g. the nondiagonal, static 
correlator S^oiq) i s of the same order as 5g (g) (22), i.e. it 
should not be neglected. Without assuming diagonality 
with respect to I, more contributions to the translational- 
rotational coupling will be taken into account which we 
expect to yield a higher value for T C MCT . If in that case 
and for l co = 2 it would be still T C MCT < T C MD , this 
would indicate that also correlators with I > 2 have to 
be taken into account, as it was done for the single dumb- 
bell in an isotropic liquid of hard spheres, however, using 
SDA |L4j . In rcf. it has been found that a one percent 
accuracy of (//™)c requires to take into account correla- 
tors up to I + 2. But on the other hand our investigation 
shows that the SDA may introduce an error of 5-10 %. 
Therefore, in order to get reliable results for the NEP, 
e.g. (/2™) c one has to remove both the diagonalization 
approximation and the restriction to I < 2. 
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